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Ultracold gases promise many applications in quantum metrology, simulation and computation. 
In this context, optimal control theory (OCT) provides a versatile framework for the efficient prepa¬ 
ration of complex quantum states. However, due to the high computational cost, OCT of ultracold 
gases has so far mostly been applied to one-dimensional (ID) problems. Here, we realize com¬ 
putationally efficient OCT of the Gross-Pitaevskii equation (GPE) to manipulate Bose-Einstein 
condensates in all three spatial dimensions. We study various realistic experimental applications 
where ID simulations can only be applied approximately or not at all. Moreover, we provide a strin¬ 
gent mathematical footing for our scheme and carefully study the creation of elementary excitations 
and their minimization using multiple control parameters. The results are directly applicable to re¬ 
cent experiments and might thus be of immediate use in the ongoing effort to employ the properties 
of the quantum world for technological applications. 


I. INTRODUCTION 

Over the last decade, the ever increasing experimen¬ 
tal toolbox of atomic, optical and molecular physics has 
lead to an exciting improvement in the control and un¬ 
derstanding of complex quantum systems [1]. Recently, 
this has resulted in an important shift of paradigm. 
While quantum systems were previously mostly studied 
to check the validity of theoretical models, interest has 
now increased in their manipulation for specific techno¬ 
logical applications. Prototypical examples for this shift 
of paradigm are atomic interferometers for quantum en¬ 
hanced metrology [2-4], atomic field probes [5] and mi¬ 
croscopes [6, 7], inertial sensors [8], atomic clocks [9], or 
applications in quantum computing [10, 11] and quantum 
simulation [12], 

In many cases, these applications rely on the controlled 
preparation of a well-defined quantum many-body state 
with particular properties. One of the key experimental 
challenges is thus the efficient transfer of a system to such 
a state. Optimal control theory (OCT) is a mathematical 
tool to devise control strategies for this transfer [13]. It 
is well studied in many physical systems, ranging from 
atoms and molecules to solid-state systems [14-18]. 

In this work, we apply it to the control of a dilute 
atomic Bose-Einstein condensate (BEC), a system which 
is well described by the three-dimensional (3D) Gross- 
Pitaevskii equation (GPE) [19, 20]. Such BECs form a 
versatile experimental platform for the storage, manip¬ 
ulation and probing of interacting quantum fields with 


’Electronic address: mennemann@acin.tuwien.ac.at 
t Electronic address: tim.langen@colorado.edu 


high precision [1], In a seminal work Hohenester et 
al. [21] demonstrated that OCT [22] provides a highly 
efficient way to realize the transfer of a BEC to a tar¬ 
get state, vastly outperforming more simple schemes. In 
this context, it has also been shown that OCT is robust 
against fluctuations and decoherence, and can also specif¬ 
ically take into account experimental constraints [23]. 
This has recently lead to first experimental demonstra¬ 
tions [24, 25]. 

OCT of BECs has so far been mostly used in one- 
dimensional (ID) settings, as the computational cost 
scales exponentially with the number of dimensions [26]. 
However, many experimental situations can only approx¬ 
imately be described by a ID model, potentially limiting 
the applicability of OCT to real-life situations. In the 
following we demonstrate the first OCT of a BEC in all 
three spatial dimensions. We go beyond situations where 
a ID approximation is feasible, thus significantly expand¬ 
ing the range of applicability of OCT for BECs. More¬ 
over, we perform an analysis of the collective excitations 
that are created as a result of the control. These excita¬ 
tions are directly connected to the non-linear nature of 
the GPE and can only be fully captured and minimized 
in a 3D treatment including multiple control parameters. 
They are thus highly relevant in realistic experimental 
situations. 


II. THE CONTROL PROBLEM 

We start with a brief review of OCT, as well as of the 
description of BECs in terms of the GPE. 
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A. Gross-Pitaevskii equation 

The mean-field dynamics of a BEC is described by the 
GPE 

K 2 

ihdti> = - — Aip+ V\ip + g\ip\ 2 ip ( 1 ) 

2m 

where ip = ip(r, 0 denotes a complex-valued wave func¬ 
tion, with initial condition ip(r, 0) = ipo £ L 2 (R 3 ;(D). 
Here, V\ = V(r,X(t)) is an external potential that 
is characterized by a single or several control parame¬ 
ters denoted by the vector A. Assuming that the wave 
function is normalized to unity, the coupling constant 
g = N4Trh 2 a s /m is defined by the mass ?n, the s-wave 
scattering length a s and the number N of atoms in the 
BEC. For example, for ultracold gases of 87 Rb atoms, 
the atomic mass is given by m = 1.44 x 10 ~ 25 kg and 
the s-wave scattering length by a 3 = 5.24 nm. Measuring 
length in units of Iq = 1 x /mi. mass in units of the atomic 
mass and time in units of to = ml^/h, equation ( 1 ) can 
be written as 

id t ip =-^Aip+V\ip + g\ip\ 2 ip ( 2 ) 

which is the starting point for our considerations below. 

B. Optimal control problem 

We seek to find an optimal time-evolution of the m- 
component control parameter 

A : (0, T) R m , A(0) = A 0 , A(T) = A T , 

which steers the system from the initial state ipo at time 
zero to a desired state ipd at final time T. Without loss 
of generality we assume that ipo and ipd are ground state 
solutions of the stationary GPE corresponding to the 
smooth external potentials V\ 0 and V\ T at times t = 0 
and t = T, respectively, with fixed parameters Ao, A t- 
To find the time evolution we apply well-known tech¬ 
niques from optimal control theory [ 21 ]. As cost func¬ 
tional, we use 

JM) = \ (1- \W d MT))\ 2 ) + l£\d t m 2 dt, (3) 

where ( u,v ) = f R3 u(r)*v(r) dr denotes the standard 
scalar product of u,v G L 2 ( R 3 ;(D). The definition 
(3) is the generalization of the functional used in Refs. 
[21, 23, 26-30] to a multi-component control parameter 
A. The first term in J measures the proximity of ip to 
the desired state ipd at the end of the steering process. 
The expression F(ip) = l — \(ipd, V ’)| 2 is known as the infi¬ 
delity and provides a measure for the difference of ip and 
ipd- In detail, it quantifies the T-norm of ip ’s component 
that is orthogonal to ipd■ The second term regularizes the 
control trajectory to account for the fact that parameters 


can never be changed infinitely fast in a real experiment. 
Here, 7 > 0 sets the penalty for fast variations of A(f). 
For our examples below we find that already a very small 
value 7 = 1 x 10~ 6 yields a satisfactory regularization. 

Our goal is to minimize J(X,ip) subject to the con¬ 
straint that ip solves the GPE (Eq. (2)) with the initial 
condition given by the respective ■ ipo . To this end, one 
introduces the Lagrange function 

L(\,ip,p) = J(X,iP)+ (4) 

Re I j p* (id t ip + ^A'0 — V\ip — g\ip\ 2 ip) dtdr 
Jo J R3 2 

where p{r,t) acts as a generalized Lagrange mul¬ 
tiplier [ 22 ]. At a local minimum (A ,ip,p) of 
J, all three variational derivatives D p L(X 1 ip,p)[Sp\ 1 
D^,L(X,ip,p)[Sip] and D\L(X, ip,p)[SX] vanish for all ad¬ 
missible variations Sp, Sip and (5A, respectively. The 
corresponding three conditions constitute the optimality 
system 


id t ip = 

-^Aip + V x ip + g\ip\ 2 ip, 

(5a) 

id t p = 

-\^P + V *P + 2 9\^\ 2 P + S’V’V, 

(5b) 

II 

CM 

— Re (ip, {d x V x )p), 

(5c) 

together with 

the initial and terminal conditions 



O 

-Se¬ 

ll 

S[ 

is 

( 6 a) 


ip(T) = -( ip d ,ip{T))ip d , 

( 6 b) 


A(0) = A 0 , A (T) = A t . 

( 6 c) 


In general, no analytical solutions are available for (5) 
with ( 6 ). Here we use an iterative method to find a nu¬ 
merical approximation of the solution. For this purpose 
it is useful to introduce the reduced cost functional 

J(X) = J(X,ip x ), (7) 

where ip\ denotes the unique solution of the Gross- 
Pitaevskii equation for a given control parameter curve 
A. The goal is to find a local (or, preferably, even global) 
minimizer A* of J. 

The most straight-forward iterative procedure that can 
be employed is the method of steepest descent, 

X k+1 =X k -a k VJ{X k ), k = 0,1,2,.... ( 8 ) 

To determine an appropriate step size a k , we perform a 
line search in each iteration: 

a k = argrnin J(A fc — aVJ( X k )). (9) 

a 

Here the upper index denotes the iteration step. A com¬ 
ment is due on the use of the gradient VJ(A fc ) in ( 8 ). 
Recall that the gradient of J at A with respect to a 
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specific inner product {-,-)x on the space X of admis¬ 
sible variations <5A is the uniquely determined element 
VJ € X such that ( X7J,SX)x = -D*J(A)[<$A] for all ad¬ 
missible variations S A £ X. The gradient thus depends 
sensitively on the choice of the inner product (-,-)^ on 
X. It has been pointed out already in Ref. [27] that any 
admissible variation (5A must have a finite value in the 
penalty term, i.e., its weak time derivative dt.5X must be 
square-integrable on (0,T), and must respect the bound¬ 
ary conditions in ( 6 c), i.e., dA(0) = SX{T) = 0. A natural 
choice for (•,-) x is thus the Hq (0, T, R m )-scalar product, 


(u,v) := f d t u(t) ■ d t v(t) dt. ( 10 ) 

Jo 

A calculation, which we present in the appendix, shows 
that this choice of (•, -)x yields 

[V J(A)] = 7 A + (d x V x )p), (11a) 

[VJ(A)](0)=0, (lib) 

[VJ(A)](T)=0, ( 11 c) 


wherein if) and p are solutions of (5a) and ( 6 a) or (5b) 
and ( 6 b), respectively. By definition, VJ vanishes at the 
boundaries t = 0 and t = T, and so the iteration ( 8 ) pre¬ 
serves the boundary conditions ( 6 c). We emphasize that 
the seemingly canonical choice of (-,-)-Y as the standard 
L 2 -scalar product would not allow to specify boundary 
data for VJ, which would result in a severe loss of sta¬ 
bility of the optimization algorithm. 

C. Implementation 

In the situations considered below we found that the 
method of steepest descent (see Eq. ( 8 )) works reliably. 
However, using more advanced methods the number of 
iterations needed to ensure convergence of the algorithm 
can be reduced significantly. In fact, our solver is based 
on the non-linear conjugate gradient scheme of Hager and 
Zhang [31], which has also been employed in Ref. [27] for 
optimal control of the one-dimensional GPE. We stress 
that all inner products and norms related to the non¬ 
linear conjugate gradient scheme need to be expressed in 
terms of the inner product given in Eq. (10). 

The reduced cost functional (7) needs to be evaluated 
several times per iteration. Moreover, at the beginning of 
each iteration a gradient vector needs to be determined 
using Eq. (11). Solutions to the time-dependent GPE 
(5a) and the adjoint equation (5b) are obtained via the 
time-splitting spectral method [32]. Initial and desired 
final states for a given potential are found by imaginary 
time propagation. 

In order to accelerate the solving of the optimal con¬ 
trol problem we perform all computations on the graphics 
processing unit (GPU) of a powerful graphics card. To 
speed up the calculations and ensure convergence of the 


algorithm we start each optimization with a coarse spa¬ 
tial grid and a relatively big time step At. The result 
for A is used as an input for another round of optimiza¬ 
tion on a finer grid. This procedure is repeated until 
the algorithm converges to a final time-evolution for A. 
A detailed description of our implementation is given in 
Appendix B. 

III. EXAMPLES 

In the following we demonstrate the results of our 
scheme by considering three applications of increasing 
complexity, which are directly connected to recent ex¬ 
periments. 

A. Harmonic oscillator potential 

In the first application we study a Bose-Einstein con¬ 
densate in an elongated harmonic potential. Initially, 
the trap frequencies are chosen such that the condensate 
is aligned along the y direction. Using a suitable time- 
evolution of the trap frequencies, we aim to rotate the 
condensate by 7r/2, while keeping it in the ground state 
of the external potential. 

An example of the transition is visualized in Fig. 1. It 
can be understood as a toy example of a broad class of 
experimental protocols in which the trapping geometry 
is changed, e.g. to mode match different traps [34], to 
(de)compress a trap [35] or to transfer condensates into 
dynamical potentials for atomtronics [36, 37]. Conceptu¬ 
ally similar pulsed manipulations are also performed to 
focus BECs in time-of-flight expansion [38]. 

1. Trapping potential 

The harmonic potential in this example is given by 

V\{x,y,z) = y (]w x (Ai)] 2 x 2 + [ijj y (X 2 )] 2 y 2 +ut 2 z z 2 ^, 

wherein the frequencies lu x and oj y can be set indepen¬ 
dently via the control parameters Ai and A 2 . More pre¬ 
cisely, we transform the external potential from an initial 
configuration with u) x = w* and u> y = w* at time t = 0 to 
a final configuration with w x = and u y = tv? at the 
final time t = T. To this end, we parametrize u) x and u> y 
as 

^x(Ai) oj x T Ai (ca x oj x ), 

A2) OJy T A2 {oJ y 

with 

Ar(0)=0, Ar(T) = 1, 

A 2 (0)=0, A 2 (T) = 1 . 
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t = 0 ms 


t = 3 ms 




t = 6 ms 


t = 9 ms 


FIG. 1: Two-parameter optimal control of an elongated harmonic potential. The timescale of the control is T = 9 ms. Initially 
aligned along the y-direction, the condensate is dynamically transformed to be aligned along the ^-direction. The black 
isosurface corresponds to the external trapping potential that is controlled using OCT, the blue isosurfaces visualize the atomic 
density. Note that for clarity only the lower half of the potential is shown. Also, here and throughout this work any trivial 
potential offset has been removed for simplicity and easier visualization. Its only effect is an overall phase shift of the wave 
function which is of no relevance to the optimization procedure. Animations of the full dynamics are available online [33]. 


We note that these parametrizations, as all others dis¬ 
cussed below, are chosen as an example and can easily 
be adjusted to the parameters accessible in a specific ex¬ 
perimental realization. 

2. Numerical simulations 

In the following simulations the number of atoms is 
N = 5000, the final time is set to T = 9 ms and u) z = 
5 kHz. The initial configuration of the trapping potential 
is given by oj 1 x = 5 kHz and w* = 0.75 kHz, the final 
configuration by = 0.75 kHz and ujf = 5 kHz. 

Before we discuss the result of the optimal control al¬ 
gorithm we first consider a numerical simulation as a 
benchmark, in which the control parameters Ai and A 2 
are varied linearly. The corresponding time-evolution of 
the trap frequencies co x and ui y is depicted in Fig. 2a. 

In order to investigate the overlap of ip with ipd beyond 
the end of the control we continue the time-evolution with 
A (t) = A(T) for t > T. We proceed analogously in the 
other examples. As can be seen from Fig. 2b the infidelity 
decreases only slightly until t = T and shows a strong 
oscillation for t > T. This behavior of the infidelity in¬ 
dicates that the final state differs significantly from the 


desired state ipd■ This is also strikingly visualized by ex¬ 
ample snapshots of the density at time t* = 22 ms in 
Figs. 2c-e. 


Next, we consider the result of the optimal con¬ 
trol algorithm. Using A °(t) = [0.25 sin(7rf/T) + 

t/T, — 0.25sin(7rf/T) +t/T] for t € [0, T] as a starting 
point, the algorithm converges to a solution that reduces 
the cost functional by four orders of magnitude. The 
time-evolution of the frequencies lu x and u> y is shown 
in Fig. 2f, the time-evolution of the corresponding infi¬ 
delity in Fig. 2g. It can clearly be seen that the infidelity 
strongly decreases until the end of the control at t = T. 
Moreover, the infidelity remains on a very low level for 
t > T, indicating that the desired final state has been 
reached with high precision. Consequently, the devia¬ 
tions of the density to the density of the desired state 
at time t = t* are very close to zero as can be seen from 
Figs. 2h-j. We note at this point that the evolution of the 
3D wave functions can naturally only be described here 
in limited detail. A supplementary video that visualizes 
these dynamics in greater detail is available online [33]. 
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FIG. 2: Two-parameter control of an elongated harmonic oscillator potential. The computational domain is chosen as 
([—10,10] x [-10,-1-10] x [—2.5,2.5])/rm 3 . In the finest discretization level we use 128 x 128 x 32 grid points and a time 
step size of At = 0.001 ms. Left column: without optimal control. Right column: optimal control. For details see text. 
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Sb 



Sb 


t = 0 ms 


t = 3 ms 




cP 

t = 9 ms 


FIG. 3: Loading of a toroidal trap using two-parameter optimal control. Animations of the full dynamics are available online 
[33]. 



FIG. 4: Saturation function used in the toroidal trap and 
splitting examples. 


B. Loading of a toroidal trap 


In the second application we consider the loading of a 
toroidal trap as shown in Fig. 3. Such toroidal traps have 
recently been employed to realize atomic analogues of 
electrical circuits to study superflow and hysteresis [39- 
43], 


1. Trapping potential 

The trapping potential is given by a slightly elongated 
harmonic potential and a Gaussian function centered at 
the origin of our coordinate system [44] 

V\(x, y,z) = ™( [^(Ai )] 2 x 2 + u] 2 v y 2 + ulz 2 ') 

+ Vb(A 2 ) exp(—2(x 2 + y 2 )/wl). 

In an experiment this Gaussian function could for exam¬ 
ple correspond to a red-detuned laser beam realizing a 
repulsive dipole potential. 

As illustrated in Fig. 3 we consider the transformation 
of the potential from an initial harmonic configuration 
with tu x = co x and Vo = 0 at time t = 0 to a toroidal 
configuration with oj x = u}£ and Vo = Vq at the final 
time t = T. Hence, a suitable parameterization of u> x 
and Vo is given by 

^x(Ai) = uj 1 x + — o4)> (12a) 

Vo(\ 2 ) = V 0 * X (\2), (12b) 

where 

Ai(0) = 0, Ai(T) = 1, 

A 2 (0)=0, A 2 (T) = 1. 

In Eq. (12b), \ pl a Y s the r °l e of a saturation function. 
The use of the saturation function ensures that Vo re¬ 
mains positive - and thus experimentally realizable - for 
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any possible choice of A 2 . This does not restrict the orig¬ 
inal control problem, as every experimentally realizable 
trajectory Vq( t) > 0 can be parametrized through a suit¬ 
able A 2 (t) in V^*x(A 2 (t)). In fact, we choose all param¬ 
eters for the external potential to be close to previous 
experimental realizations. However, our approach also 
allows us to optimize more general situations where the 
parametrization of the trapping potential is more com¬ 
plicated [45]. 

Similar saturation functions are commonly used in con¬ 
trol theory to realize limits on control parameters. In our 
particular case x is implemented using a piecewise cu¬ 
bic hermite interpolating polynomial (PCHIP). Its func¬ 
tional form is shown in Fig. 4. The interpolating points 
are chosen such that % always remains positive. More¬ 
over, y(0) = 0 and y(l) = 1. 


C. Splitting 

In terms of technological applications, a particular 
noteworthy realization of BECs is achieved using atom 
chips [47, 48]. On these chips micro-fabricated wires al¬ 
low the precise manipulation of BECs using static, radio 
and microwave fields. As a third application we thus con¬ 
sider the splitting of a single condensate into two iden¬ 
tical halves using such an atom chip [46]. A visualiza¬ 
tion is presented in Fig. 6. This splitting protocol has 
recently been used to study the non-equilibrium dynam¬ 
ics of ID Bose gases, revealing subtle effects, such as 
prethermalization [49-52], generalized statistical ensem¬ 
bles [53] and the light-cone-like emergence of thermal cor¬ 
relations [54, 55]. Moreover, it forms the basic building 
block for integrated matter-wave interferometers [56, 57]. 


2. Numerical simulations 


The following simulations are carried out using Vq = 
h x 30 kHz, wq = 5 /im , T = 9 ms and N = 5000. The 
frequencies w y = 2.5 Hz and lo z = 5 kHz are kept constant 
during the simulation. The initial configuration of the 
confinement potential is characterized by w* = 1 kHz and 
Vo (t = 0 )/h = 0 kHz, whereas the final configuration is 
given by u)f = 2.5 kHz and V 0 (T)/h = Vo* • 

As in the previous example we consider first the case 
where the parameters u) x and Vo are changed linearly (see 
Fig. 5a). Fig. 5b reveals that the associated infidelity 
does not drop at all until t = T. For t > T we observe 
a slight decrease of the infidelity. This can be attributed 
to the fact that, as time evolves, the density of the con¬ 
densate becomes more evenly distributed in the toroidal 
trapping potential, bringing its wavefunction closer to tpd- 
However, as can be seen from Figs. 5c-e, the final wave 
function still differs strongly from the wave function of 
the desired state after t* = 22 ms. 

Let us now discuss the result of the optimal control 
algorithm. An optimal time-evolution of the control pa¬ 
rameters is given in Fig. 5f. Intuitively this control can 
be understood as the result of two separate time-scales. 
During the first halve of the control, the trap frequency 
w x is increased, while the limits imposed on A 2 prohibit 
any change of Vo. During the second halve, on the other 
hand, Vo is adjusted to its final value, while lo x is only 
subject to small corrections. 

Until the end of the control this leads to a drop in the 
infidelity by approximately three orders of magnitude, as 
visualized in Fig. 5g. Furthermore, the infidelity remains 
bounded by 3 x 1CD 3 for t > T, which is well below the 
measurement sensitivity in typical experiments. Conse¬ 
quently, only slight deviations from the desired wavefunc¬ 
tion at time t* = 22 ms can be observed in Figs. 5h-j. 


1. Trapping potential 


In the experiments the splitting is realized by dress¬ 
ing the static magnetic trapping potential with a strong 
near-field radio-frequency (RF) field. The unsealed static 
potential is given by Viatic = gFTB'BiF |B|, with the mag¬ 
netic field B = (B x , B y , B z ) being well approximated by 
the famous Ioffe-Pritchard form 


B-2 

B x = B-iX - xy 
B 2 

B z = -BiZ - —zy 


By ~ Bq + 


B 2 

2 



The parameters are given by Bq = huio/rriFgFTB, B\ = 
\Jmijj\BoIrriFgpT b and B 2 = rmo^/niFgFTB- In the 
following simulations we consider s 'Rb atoms which are 
trapped in the 5 S 1 / 2 F = 2 ,mp = 2 state where gp = 
1/2. The trap parameters are 


u>o = 2-7T x 390 kHz, 
co x = cu z = ui _l = 2-7T x 2 kHz, 
u) y = W|| = 2-7T x 85 Hz. 

The resulting dressed-state potential is given by [58] 


V x = g F yBrn. F 



IB! 


+ 


B 


RF .L 


(t) 


2 


= gF/J-BrriF 



with rriF = 2, lorf the frequency of the RF radiation and 
Brff denoting the component of the linear polarized 
dressing field B^ that is aligned perpendicular to the 
static field. As in [54] we use a detuning of Arp(0) = 
—27r x 30 kHz from the mj? = 2 —»• itif = 1 transition for 
the simulation. The Rabi-frequency is parameterized by 
the control parameter A 


fWA) = u: abi X(A), D r * abi = 2^ x 155 kHz 
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FIG. 5: Loading of a toroidal trap using two-parameter control. The computational domain is given by ([—8,8] x [—8,+8] x 
[—2.5, 2.5]) ym 3 . In the finest discretization level we use 128 X 128 X 40 grid points and a time step size of At = 0.001 ms. Left 
column: linear variation of the control parameters. Right column: optimal control of the control parameters. For details see 
text. 
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FIG. 6: The splitting of a Bose-Einstein condensate, as realized by a radial deformation of an initially harmonic potential into 
a double well [46]. The two gases in the final picture are completely decoupled, with no more overlap between the respective 
wave functions. Animations of the full dynamics are available online [33]. 


wherein 

A 1 (0)=0, Ai(T) = 1. 

The control parameter A mimics the situation in exper¬ 
iments, where the double well potential is controlled by 
changing the RF field amplitude through an RF current 
in a wire. For A = 0 we recover the static harmonic po¬ 
tential, whereas A = 1 corresponds to a fully separated 
double well with no wave function overlap between the 
two halves of the system. Since the Rabi-frequency is 
strictly positive in experiments we employ the same sat¬ 
uration function \ as in the previous example (cf. Fig. 4). 

As the trapping potential is significantly changed dur¬ 
ing the splitting the atoms are radially displaced from 
their equilibrium position in the harmonic trap. Conse¬ 
quently, strong dipole and breathing oscillations are usu¬ 
ally observed in experiments. This poses a strong limi¬ 
tation to the use of such systems as interferometers [56]. 
The minimization of such excitations is therefore one of 
the main motivations for our optimization. 

2. Numerical simulations: single-parameter control 

We illustrate the splitting procedure for N = 2000 
atoms and T = 6 ms. 


In a first step we again consider the case where the 
Rabi-frequency is increased linearly (see Fig. 7a). This 
procedure is identical to the one that is typically used in 
experiments [49, 53]. At the final time t = T the infidelity 
has only decreased slightly as can be seen from Fig. 7b. 
Moreover, the infidelity shows the expected strong oscil¬ 
lations for t > T. A snapshot of the density at time 
t* = 22.5 ms is illustrated in Figs. 7c-e, revealing that 
there is large discrepancy between the computed state tp 
and the desired state ipd- 

Next, we consider the result of the optimal control al¬ 
gorithm. We find that, irrespective of the specific choice 
of A 0 , the algorithm always converges to approximately 
the same minimizer of the cost functional. The corre¬ 
sponding time-evolution of the Rabi-frequency is shown 
in Fig. 7f. We observe that the Rabi-frequency remains 
zero for the first few milliseconds. In fact, only about 
three milliseconds of the optimization time T are used 
for the transformation of the external potential. This be¬ 
havior persists even if we increase the optimization time 
T, with the Rabi-frequency vanishing for an even longer 
initial period of time. The precise timescale depends on 
the parameters of the trap, as the optimization algorithm 
tries to find a compromise between longitudinal and ra¬ 
dial directions. 


10 


f) 



d) 





- 20 ' 


30 




xi 

o' 

CO 

-ii 



10 u 


10 “ 


10 “ 


i) 


“i—i—■-1-■-r 


i - \{^dA{t ))\ 2 

\{ip d ,Si/j 2 {t))\ 2 




30 


- 20 


- 10 



N 

K 

.s 

o 

o' 


W 

.B 

o' 

S* 

co 


w 

fl 

Tj" 

o' 

CO 


FIG. 7: Splitting of a BEC using single-parameter optimal control. Left column: linear variation of the control parameter. 
Right column: optimal control of the control parameter. We note that the \('ipd, ^ 2 (t))\ 2 in (g) has been scaled and slightly 
shifted in time to account for the unknown phase and amplitude of the excitation. The computational domain is given by 
([—4,4] x [—15,+15] x [—2, 2])/xm 3 which is discretized by 96 x 128 x 48 grid points in the finest discretization level. The 
corresponding time step is At = 0.001 ms. For details see text. 
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Interestingly, our 3D control qualitatively resembles 
the result of a previous ID optimization that included 
beyond mean-field effects to model the distribution of 
atoms into the two final gases on the quantum level [57]. 
In both cases, the initial BEC is first rapidly split into 
two halves. Subsequently, these two halves are kept close 
enough to experience a tunnel coupling for a finite time- 
scale. This qualitative observation is very interesting, as 
reducing relative number fluctuations can help to signifi¬ 
cantly enhance the sensitivity of such interferometers. A 
detailed study of how useful our control can be in this 
context will be a natural extension of this work. 

As a result of the optimal control algorithm the infi¬ 
delity at the final time T is reduced by more than two 
orders of magnitude (see Fig. 7g). However, for t > T we 
again observe a strong oscillation. Snapshots of the den¬ 
sity distribution at t* = 22.5 ms are given in Figs. 7h-j. 


3. Bogoliubov-de Gennes analysis 

Interestingly, the 6 ms period of the very regular infi¬ 
delity oscillation shown in Fig. 7g for the optimized split¬ 
ting is approximately the same as the period of the infi¬ 
delity oscillation depicted in Fig. 7b for the simple linear 
splitting. This suggests that the character of the oscilla¬ 
tion is determined by the intrinsic properties of the BEC 
rather than by the splitting protocol. 

Indeed, we demonstrate in the following that the oscil¬ 
lations are caused by collective excitations of the BEC, 
which are created during, but irrespective of the details 
of the splitting process. To this end, we show that they 
are the result of a small deviation Sip from the desired 
state ipa, which can be described within the Bogoliubov- 
de Gennes (BdG) framework. 

Let therefore <I>(r,t) = cp(r)e~ ltLt l h denote an eigen¬ 
state solution of the CPE. Here, g is the corresponding 
chemical potential and <p is a solution of the stationary 
GPE, Hq4> + g\<P\ 2( P = g4>i with Hq = —H 2 /2mA + V. 
We consider a generic state ip which deviates from the 
eigenstate solution by a small fluctuation Sip, i.e, 

ip(r,t) s=s $(r,£) + 5ip(r,t). (13) 

In a linear approximation (with respect to Sip ) this small 
deviation is given by 

Sip(r,t) = ( u{r)e~ iuit + n*(r)e iu ’* t )e- i/it / ft (14) 

where u, v and u> are defined via the solutions of the BdG 
equations [59, 60] 


Ho - H ± 2 g\(p\ 2 

9<p 2 

u 


u 

-{g<P 2 Y 

-H 0 + n - 2g\cp\ 2 _ 

V 

— rujj 

V 

( 


We want to investigate small fluctuations Sip corre¬ 
sponding to some of the lowest energy eigenvalues hw in 
equation (15). To this end, we proceed in a conceptually 
similar way to [61] where numerical methods are used to 


investigate the stability and decay rates of non-isotropic 
attractive Bose-Einstein condensates. Like in Ref. [61] 
we consider the full three-dimensional problem. However, 
for the discretization of the operators in (15) we employ 
a higlr-order finite difference discretization rather than 
working in a Fourier basis. By gradually increasing the 
spatial resolution of the finite difference discretization we 
are able to verify the convergence of the algorithm. A de¬ 
tailed description of our implementation is again given in 
the Appendix. 

As an example, we find that the first three eigenvalues 
converge towards uq = ±314.54 Hz, uq = ±523.49 Hz 
and uq = ±734.26 Hz. Subsequently, the corresponding 
eigenfunctions () are normalized according to the 
norm [60] 


/ {^(r)-v 2 (r))dr= 1. (16) 

J R 3 

Knowing the frequencies uq and amplitude functions 
Ui and Vi, it is possible to investigate the time-evolution 
of the excitations given by Eq. (14). It turns out that 
\Sipi{t)\ 2 can be well described by a simple periodic os¬ 
cillation in amplitude, while the shape remains mostly 
unchanged (see left column in Fig. 8). As Ui and v, 
are purely real-valued functions, which approximately 
fulfill Vi = —Ui (see right column in Fig. 8) we find 
Sipi(r,%/2) ss Sipi(r,Ti) and hence the effective oscil¬ 
lation periods are halved with respect to the eigenvalues 
found above, i.e. T e s,i = 7)/2 = it /uq. In detail we find 
%s,i = 9.99 ms, T eff ,2 = 6.00 ms and T e s ,3 = 4.28 ms. 

Note that the effective period of the second excita¬ 
tion is very close to the period of the oscillation of the 
infidelity observed above. Indeed, plotting the time- 
evolution of \{ipd, Sip 2 (t))\ 2 along with the time-evolution 
of the infidelity in Fig. 7g demonstrates clearly that the 
oscillation of the infidelity is dominated by the second 
excitation. As further evidence, we extract the deviation 
of ip from ipd from our simulation. A comparison shows 
again very good agreement with the time-evolution of 
Sip 2 (t) (see Appendix). 

The fact that only the second but not the first excita¬ 
tion contributes to the observations can be understood 
from symmetry arguments. The first excitation corre¬ 
sponds to an antisymmetric wave function with respect 
to the longitudinal direction, whereas the second exci¬ 
tation is symmetric. During the splitting process, the 
halving of the atom number in each of the two gases, 
as well as an overall change in the longitudinal trapping 
potential leads to a symmetric change in the extension 
of the BEC in this direction. If the control is unable 
to compensate for this change in extension, the second 
Bogoliubov-de Gennes mode is automatically excited. 

This effect is especially pronounced for the linear split¬ 
ting. In contrast to that, the optimal control algorithm 
can still reduce the infidelity at t = T , but even a small 
deviation of the wave function from the stationary state 
leads to a strong oscillation in the infidelity for t >T. 
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FIG. 8: Solutions of the Bogoliubov-de Gennes equations using a 6th-order finite difference discretization for N = 2000 atoms. 
Left: Density of the first three (scaled) excitations 5i/ji(r, t), 5ip2{r, t) and 5ip3(r, t) at t = T c s,i/2, t = T a s,i/2 and t = T e «,3/2. 
Right: Normalized (with respect to the inner product (16)) amplitude functions u and v evaluated along the longitudinal 
direction at x = x s and z = 0. All functions are purely real-valued. 


Once the wave function differs from the stationary 
state in the longitudinal direction it is impossible to stop 
the observed oscillation by a simple variation of the Rabi- 
frequency. The BEC will thus oscillate for t > T after 
the end of the control. 

A central role in this scenario is played by the longitu¬ 
dinal frequency ui y . The smaller oj v the longer the exten¬ 
sion of the condensate in the longitudinal direction. In 
analogy to a classical harmonic oscillator this increases 
the susceptibility to small deviations from the equilib¬ 
rium position. We have confirmed this intuition with 
additional simulations, finding an even more pronounced 
excitation of the second mode for smaller u> y . 


This is particularly noteworthy with respect to experi¬ 
ments studying BECs in the one-dimensional limit, where 
a j y [53]. Intuitively, such experiments should be 
very well described through a ID approximation, where 
only a reduced GPE for the x-direction has to be consid¬ 
ered (see Appendix). Our results here show that such an 
approach will, in general, also lead to a strong breathing 
oscillation. Even if the ID control is able to reach the ID 
desired state with high precision, it does not necessarily 
describe the experimental reality and will thus fail in 3D. 
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4- Numerical simulations: two-parameter control 

In the last part of this article we will show how the os¬ 
cillations reported above can be eliminated using a more 
sophisticated control scheme that is made possible by the 
3D character of our control and that involves a manip¬ 
ulation of the trapping potential along the longitudinal 
direction. In experiments on atom chips, this manipula¬ 
tion can, for example, be realized using additional wire 
structures, which provide longitudinal confinement inde¬ 
pendent of the main radial trapping structures [62]. 

In analogy to the previous examples, we consider the 
following parameterization of D ra bi and lo v : 

firabi(Ai) = fi; abi x(A 1 ), ft r * abi = 2^x155kHz, (17a) 
uj y ( A 2 ) = w*A 2 , w* = 27tx85Hz, (17b) 

with 

Ar(0)=0, Ar(T) = 1, 

A 2 (0) = 1, A 2 (T) = 1. 

The only difference to the previous example is thus that 
the value of the longitudinal trap frequency u> y is now 
part of the control. We still fix oj y (t = 0) = 2w x 85 Hz 
and u) y (t = T) = 2w x 85 Hz such that the initial and 
desired final states remain unchanged. 

Using A°(t) = [■ t/T , 1] for t € [0,T] as an initial guess 
the optimization algorithm converges to a solution which 
reduces the cost functional by more than three orders 
of magnitude. The time-evolution of the corresponding 
physical parameters is given in Fig. 9a. As can be seen 
from Fig. 9b the infidelity remains very low for t > T. 
Snapshots of the density distribution at time t* = 22.5 ms 
confirm that the deviation from the desired state is ex¬ 
tremely small, see Figs. 9d-f. 

In the given example we have chosen T = 71 s, 2 - In 
contrast to that, for a time T < 7)>ff, 2 we find significantly 
worse results. The minimum time scale T is thus set by 
the oscillation period of the excitation that the control 
aims to stop. This oscillation period is in turn set by 
the geometry of the trap. Each different experimental 
situation will thus require carefully chosen parameters 
for the control. 


IV. CONCLUSION AND OUTLOOK 

In this work we have presented the first optimal con¬ 
trol of the GPE in 3D. As we have shown, this situation 
is inherently more difficult than the optimal control of 
the ID GPE because of the non-linear coupling of differ¬ 
ent coordinate directions. We have performed a detailed 


analysis of the resulting small excitations, which we were 
able to minimize by extending previous control schemes 
from a single to a multi-parameter control. 

In contrast to ID approximations our 3D approach al¬ 
lows the study of realistic trapping potentials, which will 
have direct impact on the quality of experiments and 
therefore provide an important step in the ongoing ef¬ 
fort to use the properties of the quantum world for real 
life applications. Importantly, our scheme is not limited 
to the examples discussed in this work but rather very 
flexible, with many more applications conceivable. 

A straight-forward extension of our numerical solver 
could include the treatment of excited states. This would 
allow the three-dimensional study of a recent experiment, 
where the BEC was transferred to the first excited state 
of the trapping potential via a ID optimal control se¬ 
quence [24]. Based on our observations we expect an 
even stronger excitation of BdG-Modes in such an exper¬ 
iment. In that context, another interesting application 
would be to replace the cigar-shaped confinement poten¬ 
tials used in the splitting and vibrational state inversion 
experiments by torus-shaped trapping potentials. Due 
to the different topology the issues related to the excita¬ 
tion of small perturbations are expected to be strongly 
reduced. 

Another obvious extension of this work could be to con¬ 
sider different cost functionals. More precisely, it would 
be interesting to investigate whether it is possible to re¬ 
duce the optimization times T by using other cost func¬ 
tionals which are not based on the infidelity but rather 
on a conserved quantity like the total energy. 

Finally, interesting further directions include the 
study of beyond mean-field effects using the rnulti- 
conhgurational time-dependent Hartree framework for 
bosons [63] or the optimization of finite temperature 
states. 
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APPENDIX 

A. Gradient of the reduced cost functional 

In the main text we have introduced the cost func¬ 
tional J(X,ip) in (3) and the reduced cost functional 
J(A) = J(A, if> x ) in (7). Recall that, for a given con¬ 
trol A : (0, T ) —> R m satisfying the boundary conditions 
(6c), ipx is the solution to the initial value problem (5a) 
and (6a) for the GPE with the corresponding potential 
V\. Below, we argue why the TT^-gradient of J is given by 
the component A : (0 ,T) —» R m of the solution (if>,p, A) 
to the system consisting of (5a), (5b) and 

^(A-yA) =Re(p, {d x V x )ip), (18) 

subject to the initial and terminal conditions (6a), (6b) 
and 

A(0) = A (T) = 0. (19) 

Before discussing the gradient, we first calculate the 
variational derivative of J. As it is customary in the 
context of optimization problems, we express the validity 
of the GPE (5a) in the form of a constraint Z = 0, with 
the contraint functional 

Z(\,ip) =id t ip+ ^ A^> — V x ip - g[i/j\ 2 ip- 


By definition, tp x satisfies Z{X,rp x ) = 0, hence 

J(X) = L(X,if, x ,p), (20) 

where L denotes the Lagrangian which was clehned in 
Eq. (4) of the main text. 

I/(A, ip,p) = J(A, ip) + Re f (p, Z(X,ip)) dt. 

Jo 

Eq. (20) holds for arbitrary smooth functions p : 
(0 ,T) — > L 2 (R 3 ;(D). For fixed p, differentiation of J 
in the direction dA yields 

£> a J(A)[<5A] = D x L(X,ip x ,p)[6X ] + D^,L(X, ip x ,p)[5ip] 

= D x J(X,ip x )[6X] + Re [ {p,D x Z(X, i/; x )[SX}) dt 

Jo 

+ D^J(X, ip x )[5il>] + Re [ (p,D^Z(X,^ x )[5ip]) dt, 

Jo 

(21) 

where Sip is the variation in ip x induced by the vari¬ 
ation <5A of A, i.e., it satisfies D x Z(X,ip x )[6X] + 
D^Z{X 1 ip x )[5ip] = 0 and 5ip( 0) = 0. For simplification 
of D X J , we choose p, which has been arbitrary up to this 
point, such that the last two terms in (21) cancel. Indeed, 
taking p as a solution to the terminal value problem (5b) 
and (6b), it follows that 


A/,J(A, ip x )[6ip] + Re J {p,D^Z(X,ip x )[Sip]) dt 

= -Re((ip d ,ip x (T))*(ip d ,Sip(T))) + Re f ( p,id t 5ip + Sip - V x 6ip - 2g\ip x \ 2 Sip - gip\8ip*) dt 

Jo * 

= Re J ( idtp + ^A p - V x p - 2g\ip x \ 2 p - gip 2 p *, Sip) dt = 0. 


To arrive at this result, we have performed an integra- the initial condition (6a). In view of these cancellations, 
tion by parts with respect to time, using the terminal equation (21) simplifies to 
condition (6b) and the fact that 5i/j(0) = 0 thanks to 


D x J(X)[6X] = 7 / d t X- dt(6 A) dt- Re / <p, d x V x ■ (6 A) if,) dt. 


( 22 ) 


We are now in the position to calculate the if ■'■-gradient 
of J. Recall that the Sobolev space I7 1 (0,T;R m ) con¬ 
sists of all square integrable functions A S L 2 (0,T;R m ) 


that possess a weak derivative d t X S L 2 (0, T ; R m ). Func¬ 
tions A € H 1 ( 0,T;R m ) are actually Holder continuous, 
and therefore, they have well-defined boundary values at 
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t = 0 and t = T. It is natural to consider the reduced 
cost functional J as defined on £7/(0,T;R m ), which is 
the affine subspace of functions A G H 1 ( 0,T;R m ) that 
satisfy the boundary conditions (6c). Indeed, any ad¬ 
missible control A : (0 ,T) -* R m must produce a fi¬ 
nite value in the penalty term in J, which implies that 
d t X G T 2 (0, T;R m ). The tangent space to £1/(0, T;R m ), 
i.e., the space of possible variations <5A, is the linear sub¬ 
space ££/(0, T; R m ) of all functions A G H 1 (0,T;W n ) 
with vanishing boundary values, A(0) = A(T) = 0. This 


is a Hilbert space with respect to the inner product 

(Ai, A 2 ) := f ■ d t A 2 (t) dt. 

Jo 

By definition, the gradient of J with respect to the inner 
product (•,•) is the uniquely determined element A G 
i£Q(0,T;R m ) such that (A, <5A) = DaJ(A)[ 6A] for all 
variations <5A G ££/(0, T; R m ). In view of (22), A satisfies 


f d t (A — yA) • d t (£A) dt = — f Re(p,d\V\ip) ■ (6A) dt for all variations <5A G £1/(0, T; R m ), 

Jo Jo 


(23) 


and A G £1/(0, T;R m ) induces the boundary condi¬ 
tions (19). To verify that the solution A to the boundary 
value problem (18) and (19) satisfies (23), it sufficies to 
integrate by parts in the time integral on the left-hand 
side, using that <5A(0) = SX(T) = 0. 

B. Algorithms and implementation 

1. Numerical evaluation of the cost functional 

The evaluatation of the reduced cost functional (7) for 
a given control curve A implicitly involves the computa¬ 
tion of ip\, that is, the solution of the GPE. No analytical 
solutions are available in general, so we use a numerical 
approximation. For brevity of notation, we write ip in¬ 
stead of ip\ in the following. 

For the numerical computation of the first term in 
(3), that is 1/2 (l — \(ip d ,ip(T))\ 2 ) , we have to solve the 
GPE (5a) with initial data (6a) for t G [0, T]. Our simu¬ 
lations are performed on the spatial domain 

II = [—L x /2, L x /2] x [—L y /2, L y /2] x [~L z /2,L z /2] 

with L x , L y and L z chosen sufficiently large to capture 
the significant part of the rapidly decaying solution ip. 

For numerical discretization in time, we employ the 
following time-splitting spectral method [32]: 

^(t n+ i) » e- iB Z At / 2 e- iAAt e- iB * At / 2 iP(t n ), (24) 

with operators A = — 1/2A, B ± = V\(n+i/2) + g\ip A \ 2 , 
and with t n = nAt, n = 0, ...,N — 1 s.t. NAt = T. 
Here X n+1 /‘ 2 '> = 1/2 (A (t n ) + A(i( n+1 ))), and the choice 
of ip// is given below. Thus, the nth time step con¬ 
sists of the following three sub-steps. First, solve idtip = 
(Vx<.n+ 1 / 2 ) +g\ip(t n )\ 2 )ip for a duration of At/2 with initial 
value ip(t n ); thus ip~ = ip{t n ). The result is used as initial 
value for the free Schrodinger equation id t ip = —1/2 A ip, 
which is then solved for duration of At; the result is ip+. 


Finally, id t ip = (P A( „ + i/ 2 ) + g\ipn\ 2 )ip is solved with ini¬ 
tial value ip+, again for a duration of At/2. The result 
of the third sub-step is taken as ip(t n + 1 ). 

The free Schrodinger equation is solved using the 
Fourier spectral method. To this end, the wave func¬ 
tion ip is interpolated by a trigonometric polynomial on 
the grid points of the cartesian grid 

( X jx 1 Vjy > Z N ) = 

{~J j x/‘J‘ + jx Ax, —Ly/2 + j y Ay, —L z /2 + j z Az), 

where Ax = L x /J x with j x = 0,..., J x — 1 etc. Thus, at 
time t n , the wave function ip(t n ) is represented by a three- 
dimensional array of complex numbers ip^ G 

As Matlab code, the ?rth time step looks as follows: 

ip = exp(— li * (V A („ + i/3) + g * abs(ip)A2) * At/2) .* ip ; 
ip = fftn(i/>); 
ip = M .* ip; 
ip = ifftn( , 0); 

ip = exp(— li * (V^C"+ 1 / 2 ) + g * abs(ip).~2) * At/2) .* ip; 

(25) 

The array M represents the action of the free 
Schrodinger operator. Due to the vectorized implemen¬ 
tation in Matlab this procedure is highly efficient. 

The method is of second order in time and of spec¬ 
tral accuracy in space, provided that ipo and V\ are suf¬ 
ficiently smooth. In comparison to a finite difference 
Crank-Nicolson scheme (see for example [21, 32]), the 
solution of a linear evolution equation is avoided, and 
less grid points J = J x J y J z are needed to achieve the 
same quality of approximation for ip. 

Typically, the numerical costs for our implementa¬ 
tion (25) are dominated by the fast Fourier transforms 
f f tn and if f tn, which are of order 0( J log J). However, 
in some simulations (Splitting), the costs for computing 
the external potential V X (n.+i/ 2 ) exceed that of the Fourier 
transforms. 
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For the numerical solution of the optimization prob¬ 
lem, on the order of 10 to 100 evaluations of the cost 
functional are needed. The respective solution of the 
time-dependent GPE is performed on the graphics pro¬ 
cessing unit (GPU) of a powerful graphics card. Thanks 
to the vectorized implementation (25), it suffices to ini¬ 
tialize the arrays xp and V once at the beginning, using 
the Matlab command gpuArray. For handling the inter¬ 
mediate results or for calling the data in the memory of 
the main processor at the end of the computation we use 
the command gather. The trap potentials need to be 
updated in each time step. However, these calculations 
can be performed in a vectorized way on the GPU as 
well. 

Finally, we compute 1/2 (l — \(xp d ,xp(T))\ 2 ) with 
xp(T) « -0 (JV) , using a quadrature formula. The integral 
^ f 0 |<9 f A(t)| 2 dt is computed by a quadrature formula as 
well, using a finite difference formula of second order for 
the approximation of the time derivative (9 t A- 

2. Numerical computation of the gradient 


A further difficulty in the numerical computation of 
the adjoint equation arises from the conjugate-complex 
quantity p* in gip 2 p*. Without going into details, we 
refer to the implementation in [27], which can be easily 
applied to the three dimensional case. As in the case of 
the Gross-Pitaevskii equation the computation of r can 
be significantly accelerated by using the graphics card. 
Still, the costs for the computation of the gradient are 
three to four times higher than for the evaluation of the 
cost functional. 


3. Computation of the initial and desired final states 

The initial and terminal states xp 0 and xp d are assumed 
to be ground state solutions of the stationary Gross- 
Pitaevskii equation. We compute them by imaginary 
time propagation [65, 66] (also known as normalized gra¬ 
dient flow): the time step At in (25) is replaced by —iAt, 
and the wave function <p is normalized after every time 
step. By using adaptive time stepping, we reach a suffi¬ 
ciently exact solution with justifiable numerical costs. 


According to (11a), the ^-gradient J(A) is obtained 4 Mer detoffa gf fhe implementation 

as solution to the second order problem 


|i[VJ(A)]=r(A), 


r( A) = qA + R e(xp, (d\V\)p) 

(26) 


subject to the boundary conditions [V j(A)] (0) = 0 and 
[V J(A)] (T) = 0 . The time derivatives are discretized by 
second order finite differences. 

To evaluate the right-hand side r(A), the functions if 
and p need to be determined for t € [0, T\ . First, the state 
equation (5a) is solved as described above. Then, the ad¬ 
joint equation (5b) is solved backwards in time, for the 
terminal condition (6b). For solution of the adjoint equa¬ 
tion, a time-splitting method is applied as well: we alter¬ 
nately solve the equations id t p = V\p + 2g\xp\ 2 p + gxp 2 p* 
and id t p = —1/2A p. The free Schrodinger equation is 
discretized by the Fourier-spectral method, and the value 
of d\V\ at time t = (n — 1/2)At is computed by means 
of the complex-step derivative approximation [64]. 

For integration of the ajoint equation on the time in¬ 
terval \{n — l)Af,nAt], an approximation of the wave 
function ip( n ~ 1 / 2 ') = 1/2 {xp^ n ~ v > +i/>(”)) is needed. Since 
it is impossible to store the arrays ip'A for every time 
step n = 0on the graphics card, the state equa¬ 
tion is simultaneously solved backwards in time as well. 
The procedure is sketched in Fig. 10: the calculation 
of involves only two instances of the wave func¬ 

tion and the “old” adjoint state — that is, xp ^, 
and p (Tl b As soon as the approximations of xpA- 1 ) and 
p(n-i) are available, also r^ n ~ 1 ^ can be computed. In this 
way it is enough to store at each time step four arrays in 
three dimensions, and the values of all available r ^ with 
n = 0,..., N (the storage space of which is neglegible). 


For the numerical solution of the considered optimal 
control problems we use a personal computer (i7-4770K 
CPU @ 3.50Ghzx8) and Matlab. The parts with the 
highest numerical costs, thus the solving of the partial 
differential equations and the computation of the external 
potentials, are performed on the graphics card (GeForce 
GTX TITAN), which accelerates the calculations signif¬ 
icantly. The evaluation of the Fourier transform, for ex¬ 
ample, on the finest space discretization can be acceler¬ 
ated by a factor 4-6. In this context, it is important to 
mention that the CPU-version of fftn in Matlab is par¬ 
allelized as well and hence uses all cores available on the 
CPU. 

In general it is useful to initially solve each optimal 
control problem with a small number of Fourier modes 
J x , J y , J z and with a relatively big time step At. Sub¬ 
sequently, the same optimal control problem is solved on 
a finer mesh grid and with smaller time step, whereby as 
initial data A 0 is used, obtained as approximated solu¬ 
tion in the computation before. We repeat this procedure 
until the computed control curve with respect to the old 
discretization does not differ from the control curve of 
the finer discretization anymore. 

We consider a sequence of discretization parameters 

(4”.-A A,(a*) ( i >) -» (4 2, A 2) . A.(A 2 >) 
(4">, 4">, 4*°. (a*) (M| ) 

with Jx +lS> > Jx\ > Jy C \ </i £+1 ' ) > Jz ^ and 

(At)R +1 ) < (Ai)W for t = 2,..., M. Typically on the 
order of 10 to 100 iterations of the conjugate gradient 
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FIG. 10: Computation of the source term r needed to determine the gradient: the calculation of the array p*™ -1 ) involves only 
two instances of the wave function ip^ n \ i/A" -1 ' and the current adjoint state p^K As soon as the approximations of 
and p (n_1 ^ are available, also rA n_1 ) can be computed. At each time step only the grey shaded objects need to be stored in 
the memory of the graphics card. The storage space for r il) with £ = 0 ,..., N is negligibly small. 


method are needed for solving the optimal control prob¬ 
lems on the coarse grid. The computational time is of 
some minutes. The numerical costs for the calculations 
of each single iteration increase rapidly with each dis¬ 
cretization level. In the same time if one gets near to the 
local minimum, less iterations for finding the local min¬ 
imum are required. By means of the described strategy 
each of the presented optimal control problems can be 
solved in several hours computing time with respect to 
the finest discretization level. 


C. Numerical solution of the 3D Bogoliubov-de 
Gennes equations 

For numerical treatment of (15), we proceed analo¬ 
gously to [67]: a change of variables u = \(w\ — w 2 ) and 
v = ^{wi + w 2 ) transforms the system into: 

_ H 0 - p + gc/> 2 

H 0 - p + 3gcj) 2 

(27) 

A double application of the operator decouples the eigen¬ 
value problem, 

(H 0 - (i + g<t> 2 )(H 0 - P + 3# 2 )w;i = Atui, (28a) 

(H 0 - p + 3g</> 2 )(-ff 0 - p + g(f 1 )w 2 = A w 2 , (28b) 

where A = h 2 u> 2 . Clearly, it suffices to solve the first 
eigenvalue problem (28a). 

The eigenvalue problem Aw\ = Awi given in (28a) can 
only be solved using numerical methods. To this end, 
the operator A = (H 0 — p + g(/) 2 )(Ho — p + 3 gcj) 2 ) is dis¬ 
cretized via a 6th-order symmetric finite difference for¬ 
mula. Clearly, (f> und p need to be determined in advance 
and with high precision. Here, we solve 

H 0 (f> + g\(j>\ 2 (j> = p</>, Ho = — h 2 /2m A + V 


w 1 
w 2 


= Huj 


w 1 
w 2 


using the same 6th-order finite difference discretization 
along with the imaginary time-stepping algorithm (see 
above). In this context, the second-order time-splitting 
method is replaced by the classical Runge-Kutta method 
of order 4. Subsequently, the chemical potential can be 
computed using the identity 

P = / (^|V</>(r)| 2 + I/(r)|<()(r)| 2 -(- 5 |<()(r)| 4 )dr. 

Jr 3 v 1 ' 

Once (f> and p have been determined we need to solve the 
discretized eigenvalue problem (28a). 

Naturally, we consider the same computational domain 
([—4,4] x [—15, +15] x [—2, 2]) pm 3 that was used in the 
splitting experiment in section III C. Like in the original 
experiment we employ J x = 96, J y = 128 and J z = 48 
grid points in the respective coordinate directions (in 
the finest discretization level). The resulting large-scale 
eigenvalue problem is then solved efficiently by means of 
an iterative algorithm. For this purpose we employ the 
Matlab function eigs which only determines the most 
relevant eigenvalues and their corresponding eigenfunc¬ 
tions: the algorithm yields the eigenvalues closest to a 
specified shift a which we set to a value slightly larger 
than zero. (We are only interested in the first few non¬ 
trivial solutions of (28a) corresponding to the eigenval¬ 
ues of smallest magnitude.) The underlying algorithm of 
eigs requires the repeated solution of the linear system 
of equations 


(A — al)x = b (29) 

for a given right hand side b. We employ the bicon¬ 
jugate gradients stabilized method (bicgstab) which is 
implemented in Matlab as well. Note that A — crl is 
badly conditioned which is why the bicgstab-routine 
needs to be called with a preconditioner M = M-\ M 2 , i.e. 
equation (29) is effectively replaced by M~ 1 (A — crl)x = 
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M~ l b. We found that the algorithm converges reason¬ 
ably fast when the factors M\ and M% are given by the 
matrices L and U obtained from a sparse incomplete 
Lf7-factorization. Such an approximate factorization of 
A — al can be computed using another Matlab function 
called ilu. For further information about the Matlab 
functions mentioned above we refer to the Matlab doc¬ 
umentation and the literature cited therein. The time 
needed to compute a few eigenvalue-eigenvector solutions 
of the Bogoliubov-de Gennes equations depends strongly 
on the number of grid points J x , J y and J z . For the num¬ 
ber of grid points reported above the whole computation 
takes on the order of five hours computing time utilizing 
the above mentioned CPU. 


D. Extracting the excitation from the 
time-evolution of the wave-function 



FIG. 11: Density of Aat t = 9 ms. 

The small perturbation which causes the oscillation of 
the infidelity in the splitting example can be extracted 
directly from the time-evolution of the wave function if). 
To this end, we assume that ip(r,t) and 0d(r) are almost 
identical for t = T, i.e., 

i P(t = T) « e i( Vd, 0 = argmin \\i/j(T) - e l9 'ip d \\. 

6 ' 


This assumption is in good agreement with our obser¬ 
vations, where the minimum value of the infidelity is 
reached at this point in time. In analogy to equation (13) 
we define the difference A0 := i[)(r,t) — $(r,t), which 
leads to the result 

A0(r, t) := 0(r, t ) - e i V«i(r)e- , ' i{ ‘- T)/s , t > T. 

Here, we have introduced an additional phase factor 
e i/xT/s orc j ei . take into account that we consider the 
time-evolution of A 0 starting at t = T. A snapshot of 
the density |A?/>(r, f)| 2 for t = 9 ms is shown in Fig. 11. 
It is quite obvious that the distribution of the density 
is very similar to the distribution of the density of the 
second excitation depicted in Fig. 8. 

E. One-dimensional approximation for the splitting 
of a BEC 

We briefly discuss the ID approximation for the split¬ 
ting example. In this case, the reduced GPE for the 
x-direction is given by 

ihdtij) = -|^<9 ;m ,V’ + U a (x,0,0)i/> + g ld \i/j\ 2 i/j, 

where the effective ID interaction strength gi d is found 
by integrating out the two transversal dimensions [68, 69] 

/ oo /*oo 

/ \4>{y,z)\ 4 dydz. (30) 

-oo J — oo 

Here, <p(y, z ) := 0(0, y, z) corresponds to the normalized 
ground-state solution of the 3D model in the (x = 0)- 
plane. 

With this approximation we find yid ~ h x 
1300.44 Hz ym for N = 2000 atoms. This value describes 
the situation along the whole x-axis and also leads to 
reasonable results away from the center of the cloud, as 
can be seen from Figs. 12a-b. We then follow the same 
procedures as in the 3D case to find an optimal control 
trajectory for the Rabi frequency. 
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FIG. 12: Initial state a) and desired state b) of the splitting example along the transversal ^-direction The blue solid lines corre¬ 
spond to the eigenstates of the ID approximation using the effective coupling constant gn. Dashed lines correspond to the eigen¬ 
states of the full 3D GPE evaluated along the ^-direction for shifted values of y and 2 . Here, <j> 3D,i(a:) = (/) 3 d(x, 7.5 pm, 1 pm), 
4 > 3D,2(2:) = <j> 3D(®, 7.5 pm, Opm), rj> 30,3(2:) = ^3D(i,0pm,lpm) and (p 3 o, 4 {x) = (/>3d(®, 0 pm, 0 pm). Each wave function has 
been normalized to unity, c) Optimal control of the Rabi-frequency corresponding to the ID approximation, d) Infidelity 
(ID-ID) corresponding to the one-dimensional model and infidelity (1D-3D) when the same trajectory of the Rabi-frequency 
is applied in a simulation using the 3D model. 














